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RESUME 

Dans  ce  rapport  on  dSveloppe  un  schema  numgrique  de  solution  de 
l'equation  unidimensionnelle  non  lingaire  de  conduction  de  la  chaleur. 
On  applique  ensuite  ce  schema  I  1' interaction  laser -metaux  et  on  montre 
l’importance  de  la  variation  des  propri€t€s  thermophysiques  sur  la 
distribution  de  la  temperature  dans  une  plaque  metal lique.  (NC) 


abstract 

This  report  establishes  a  numerical  scheme  to  solve  the  nonlinear 
one-dimensional  heat  transfer  equation.  The  scheme  is  then  applied  to 
the  case  of  laser-metals  interaction  and  used  to  show  the  importance  of 
thermophysical  properties  variation  on  the  temperature  distribution 
within  a  slab./fU) 
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1.0  INTRODUCTION 


The  development  of  the  COj  laser  with  its  high  efficiency  and 
high  output  power  has  made  possible  a  variety  of  laser-target  inter¬ 
action  experiments.  However,  as  the  cost  of  performing  exhaustive 
parametric  studies  to  establish  potential  applications  is  prohibitive, 
computer  simulations  were  developed  which  can  predict  laser  effects  with 
sufficient  accuracy  to  reduce  the  number  of  required  experiments  to  an 
acceptable  level. 

Although  Hanley  (Ref.  1)  has  developed  an  operational  finite 
difference  code  to  solve  the  three-dimensional  heat  flow  equation,  his 
solution  does  not  conserve  energy  and  leads  to  inaccurate  results.  For 
'simplicity,  in  this  report,  we  solve  the  one-dimensional  problem 
by  properly  treating  the  boundary  conditions  to  conserve  energy. 
Exploitation  of  the  Kirchhoff  transformation  of  the  temperature  scale 
allows  temperature-dependent  physical  properties  to  be  included  relati¬ 
vely  simply.  Furthermore,  by  reducing  the  truncation  error  to  terms  of 
second  order  in  time  and  fourth  order  in  space,  we  obtain  significantly 
more  accurate  results  for  a  specific  time  or,  conversely,  we  reduce  the 
computation  time  for  a  given  accuracy.  The  extension  of  this  numerical 
scheme  to  three-dimensional  problems  is  straightforward. 

Section  2  deals  with  the  mathematical  formulation  of  the  problem, 
and  Section  3  develops  the  finite-difference  approach  used  to  solve  the 
one-dimensional  heat  diffusion  equation.  In  Section  4,  we  present  some 
numerical  results  and  compare  them  with  the  few  existing  analytical 
solutions.  In  Section  5,  we  discuss  our  results  and  their  limitations. 

This  work  was  performed  at  DREV  between  May  1978  and  February 
1979  under  PCN  33B06,  Effects  of  Laser  Beams  on  Materials. 
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2.0  FORMULATION  OF  THE  PROBLEM 


The  mathematical  treatment  of  laser-target  interaction  involves 
the  solution  of  the  differential  equation  of  heat  conduction  in  bounded 
media  under  appropriate  initial  and  boundary  conditions.  For  a  sta¬ 
tionary,  homogeneous,  isotropic  solid  with  internal  heat  generation 
the  differential  equation  in  a  Cartesian  coordinate  is  (Refs.  2,  3) 


P  C  (T)  21  (Z,t)  =  — 
P  3t  3z 


K(T)  21  (z,t) 
3z 


+  g[z.t] 


[1] 


In  this  equation,  p  is  the  density,  C^,  the  specific  heat  at 
constant  pressure,  T,  the  temperature,  z,  the  spatial  coordinate,  t,  the 
time,  K,  the  thermal  conductivity  and  g[z,t],  the  internal  heat  source 
or  sink. 


This  differential  equation  of  heat  conduction  has  numerous 
solutions,  unless  a  set  of  suitable  boundary  and  initial  conditions 
are  prescribed.  In  this  study,  we  consider  linear  and  nonlinear  bound¬ 
ary  conditions  of  the  first  and  second  kind.  Mathematically,  they  have 
the  following  forms: 

T  =  To  or  0  [2a] 

and/or 

K[T]  21  =  F  [2b] 

3r,i 

where  To  is  a  constant  temperature,  3/3n^  denotes  differentiations 
along  the  outward-drawn  normal  at  the  boundary  surface,  S^,  and  F  can 
be  arbitrary  functions  of  time  and  surfaces. 

When  the  thermal  properties  of  the  solid  vary  with  temperature, 
the  partial  differential  equation  is  nonlinear.  If  it  is  assumed 
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that  p ,  0^  and  K  are  dependent  on  temperature  but  independent  of  position 
and  time,  and  that  the  heat -generation  term  g  =  g(z,t)  is  independent  of 
temperature,  then,  changing  the  dependent  variable  using  the  Kirchhoff 
transformation  (Refs.  2,4,5)  will  remove  the  thermal  conductivity  from 
the  differential  operator. 


This  transformation  is  accomplished  by  defining  a  new  dependent 
variable,  U,  as 


[3] 


where  U  =  U  (T) ,  T  =  T  (z,t)  and  K0  is  the  value  of  thermal  conductivity 
for  T  =  0.  Since  K  is  a  function  of  T  only  eq.  1  can  be  written  in  the 
following  form 


PC 


3T 

at 


2 

a  t 
~T 


3K 

Expressing  —  in  the  form 
3z 

3K  _  dK  3T 

3z  dT  3z 

and  substituting  this  into  eq.  4  gives 


[4] 


[5] 


K  ~ 

3z  dT 


[6] 
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From  eq.  3  we  find 

3U  =  du  3T  =  K_  3T 

3t  dT  3t  K0  3t 

3U  _  dU  3T  JC_  3T 

32  dT  3z  K0  32 

32U  =  3_  K_H._JL.Ui  31+  K  ~ 

3z2  3z  K0  3z  K0  3z  3z  3z2 

Substituting  eq.  7  in  eq.  6  gives 

1  3U  32U  g 
2  + 

a  3t  3z  Ko 


[7a] 

[7b] 

[7c] 

[8] 


Since  the  thermal  diffusivity  o  =  K/pCp  in  eq.  8,  is  a  function 
of  temperature,  the  equation  is  still  nonlinear,  but  simpler  in  form. 
Through  this  transformation,  the  boundary  conditions  become 

U  =  U0 

and/or 

Ko  »L  .  F 
3ni 

Even  with  these  simplifications,  the  solution  of  the  differen¬ 
tial  equation  is  very  complex  and  analytical  results  can  only  be  obtained 
for  a  restricted  number  of  particular  or  specific  cases,  which  we  will 
use  without  derivation  to  check  the  valicity  of  our  numerical  model. 
Following  Hanley  (Ref.  1) ,  we  use  a  standard  finite-difference  approxi¬ 
mation  to  solve  this  one -dimensional  heat  conduction  equation.  One  of 


[9a] 

[9b] 
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our  aims  is  to  obtain  a  good  compromise  between  computer  time  and 
accuracy  of  the  results,  particularly  those  related  to  a  comparison 
with  experimental  ones. 


3.0  FINITE-DIFFERENCE  MODEL 


3.1  Finite-Difference  Approximation 

We  use  an  explicit,  central  difference  scheme  to  solve  the 
following  differential  equation: 

3U  .  ..  a3^U  ,  ag 

—  (z.t)  =  - j  (z.t)  +  -fi¬ 
st  3z  K0 


[10] 


We  consider  only  materials  opaque  to  the  incident  laser  beam 
(i.e.  materials  such  as  metals  for  which  the  absorption  depth  is  much 
less  than  the  wavelength  of  the  laser  radiation)  and  chemically  as  well 
as  nuclearly  stable.  Under  such  circumstances,  there  is  no  internal 
heat  source  or  sink  and,  therefore,  g  (z.t)  =0. 

The  boundary-value  problem  under  consideration  is  that  of  a  slab 

of  finite  thickness  whose  front  surface  has  a  uniform  and  constant 

thermal  load  (i.e.  3U/3z  =  Cte  at  z  =  0),  whereas  the  back  surface  is 

3U 

either  insulated  (i.e.  =  0  at  z  =  L)  or  held  at  the  ambient  temper¬ 

ature  (i.e.  U  *  U0  at  z  =  L) . 

We  use  a  Taylor  series  expansion  of  the  function  as  our  basic 
concept  in  the  finite-difference  approximation  of  the  differential 
equation  and  the  related  boundary  conditions.  Since  the  procedure  is 
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quite  straightforward,  we  will  give  only  the  main  results  here.  When 
a  function  U(z,t)  and  its  derivatives  are  finite,  continuous,  and  have 
a  single  value,  this  function  can  be  expanded  in  the  form  of  the  Taylor 
series  about  the  point  z  as  ^ 

2  2  3  3 

U(z+Az,t)  =  U(z,t)  +  Az.  —  (z,t)  +  (z,t)  +  (z,t)+... 

3z  2!  3z  3!  3z 

HI] 


or 

U(z-Az,t)  =U(z,t) 


Az. 


3U 

3z 


2  2  3  3 

-  Az  3  U  ,  ..  Az  3  U  r  ... 

(z,tj  + - -  (z,t) - J  (z,t)  +... 


2!  3z 


3!  3z‘ 


[12] 


The  first-order  central -difference  approximation  is  obtained  by 
subtracting  eq.  12  from  eq.  11 


3U 

3z 


(z,t) 


z 


UCz*Az,t)  -  U(z-Az,t)  +  0(-Az2-j 
2Az 


[13] 


2 

The  term  0(Az  )  on  the  right-hand  side  indicates  that  the  error 
involved  in  cutting  off  the  infinite  series  is  of  the  order  of  (Az2) . 
Similarly,  the  addition  of  eqs.  11  and  12  gives  the  second  derivative 
of  the  function  as 


32U  f 
— j  (z,t) 

3z 


U(z+Az,t)  ♦  U(z-Az,t)  -  2U(z,t)  t  q(az2).  [14] 

Az2 


where  the  truncation  error  is  of  the  order  of 


(Az2). 


We  also  have  for  the  time  variable  the  following  finite-difference 
approximation 


& 
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2H.  (z,t)  =  ud.t.it)  -  u(,.t)  ,  0(4t) 

at  At 


[15] 


We  then  find,  for  the  finite-difference  approximation,  the 
following  relation 


U(R,£+1)  =  U(R,£)  +  — 

[Az] 


U(R+1,£)  +  U(R-1,£) 


2aAt 

[Az]2 


U(R,£)  +  0(At)  +  0(Az2) 


[16] 


where  Az  is  the  spatial  increment  and  At,  the  temporal  one;  R  refers  to 
lattice  points  along  the  z  axis  and  t,  to  the  integral  multiple  of  the 
step  At  along  the  time  axis.  The  coefficient  a  is  the  temperature- 
dependent  thermal  diffusivity  and  can  be  expressed  as 

«CT)  =  - ^ -  [17] 

P(T)  .  Cp(T) 

3.2  Boundary  Conditions 

Using  the  central -difference  approximation,  the  boundary  condi¬ 
tions  become 


1 .  Back  surface 

a)  Held  to  ambient 

U(Rmax,  l)  =  U„  [18a] 

b)  Insulated  (i.e.,  no  heat  flow) 

U (Rmax .  *1,1)=  UfRmax.  -1,  t)  [18b] 
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and 


2 .  Front  surface 

2  A  7 

U(Rmin.  -  1,1)  =  U(Rmin.  +  1,1)  +  — .  e  (T)  .  <j>  .  COS  ©[19] 

K0 

where  Rmax.  and  Rmin.  indicate  respectively  the  maximum  and  the  minimum 
values  on  the  z  axis,  K0  is  the  ambient  temperature  thermal  conductivity, 
e(T)  the  temperature -dependent  absorption  or  coupling  coefficient,  p , 
the  incident  flux  or  power  density  on  the  material  and  0,  the  angle  of 
incidence  of  the  laser  beam  with  respect  to  the  normal  to  the  front 
surface. 

3.3  Errors  Involved  in  Finite  Differences 


Since  in  the  process  of  the  numerical  solution  of  differential 
equations  the  derivatives  are  approximated  with  finite-difference 
expressions  at  each  nodal  point,  an  analysis  of  the  possible  errors 
involved  and  of  the  way  to  reduce  them  is  paramount.  There  are  two 
main  types  of  errors:  round-off  and  truncation  errors.  Furthermore, 
because  of  the  explicit  numerical  scheme  used,  a  specific  relation 
betw.een  the  spatial  and  the  temporal  variables  must  be  satisfied.  This 
relation  is  called  the  stability  criterion.  A  detailed  derivation  of 
the  results  given  below  is  beyond  the  scope  of  the  present  work. 

However,  interested  readers  can  consult,  for  more  information,  anyone 
of  the  following  publications  (Refs.  6-9). 

3.3.1  Round-Off  Errors 

Numerical  calculations  are  carried  out  only  to  a  finite  number 
of  significant  figures.  At  each  step,  the  error  introduced  by  rounding 
off  the  numerical  calculations  is  called  the  round-off  error.  In  linear 
problems,  the  effects  of  these  errors  superimpose  themselves  in  the 
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solution.  The  use  of  small  mesh  size,  although  desirable  for  a  better 
approximation  of  the  differential  equation,  increases  the  cumulative 
effect  of  round-off  errors.  Therefore,  one  cannot  always  say  that  de¬ 
creasing  the  mesh  size  necessarily  increases  the  accuracy  in  finite- 
difference  calculations.  On  the  other  hand,  carrying  out  the  numerical 
calculations  at  the  intermediate  stages  to  two  or  more  additional  signi¬ 
ficant  figures  will  reduce  the  cumulative  effect  of  round-off  errors  at 
the  expense  of  increasing  the  computation  time.  However,  since  the 
distribution  of  these  errors  has  many  of  the  features  of  a  random  process, 
it  is  likely  that  the  effects  of  these  errors  will  generally  cancel  each 
other  in  part.  Therefore,  it  is  very  difficult  to  determine  exactly 
the  order  of  magnitude  of  the  cumulative  departure  of  the  solution  due 
to  round-off  errors. 

3.3.2  Stability  of  Finite-Difference  Solutions 

At  each  stage  of  the  calculations,  some  round-off  errors  will  be 
present.  The  solution  of  finite-difference  equations  is  called  stable 
if  the  difference  between  the  exact  and  the  numerical  solutions  tends 
to  zero  as  At  and  Az  tend  to  zeTO  and  does  not  increase  exponentially. 
Specifically,  the  solution  will  be  stable  if  the  following  condition  is 
satisfied: 


It  should  be  noted  that  the  form  of  the  difference  equations 
depends  on  the  type  of  the  differentiation  scheme  used,  that  of 
differential  equations  and  the  boundary  conditions.  Therefore,  the 
stability  criteria  given  above  cannot  be  generalized  for  other  systems. 
In  fact,  each  system  must  be  examined  individually  for  stability. 
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Unfortunately,  there  is  no  general  method,  for  nonlinear  problems,  that 
can  be  used  effectively  to  determine  the  stability  of  the  resulting 
finite-difference  equation. 

3.3.3  Truncation  Error 

The  Taylor  series  expansion,  used  in  expressing  a  partial  differ¬ 
ential  equation  in  finite  differences,  is  truncated  after  a  prescribed 
number  of  terms.  The  error  involved  in  each  step  of  calculation  result¬ 
ing  from  the  truncation  of  the  series  is  called  the  truncation  error. 

2 

In  our  case,  that  error  involves  terms  of  order  At  and  (Az)  .  As  the 

mesh  size  is  reduced,  and  accordingly  the  time  step  to  satisfy  the 

stability  criteria,  the  truncation  error  is  expected  to  become  smaller 

so  that  the  numerical  solution  converges  faster  to  the  true  solution. 

Of  course,  this  increases  the  number  of  nodal  points  and  the  computation 

time.  It  is  interesting  to  note  that  the  truncation  error  is  reduced 

2  4 

to  terms  of  order  (At)  and  (Az)  by  satisfying  the  following  relation: 

-  [21] 

6 


a  .At 


(Az)‘ 


Under  this  condition,  the  finite-difference  solution  approaches, 
within  a  determined  residual  error,  the  true  solution  of  the  differential 
equation  at  a  faster  rate  and  it  will  be  used  in  our  numerical  scheme. 
This  value  is  called  the  "stability  constant"  and  it  will  be  referred  to 
under  this  name  later  on. 

3.4  Program  Summary 

A  computer  program  has  been  written  in  FORTRAN  IV  to  evaluate 
the  finite-difference  approximation  described  in  the  previous  section. 

A  listing  appears  in  Appendix  A  along  with  detailed  running  instructions. 
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The  program  can  deal  with  different  types  of  opaque  materials  as 
long  as  their  thermophysical  properties  and  their  temperature  variation 
are  known.  The  number  of  points  in  the  lattice  along  the  z  axis  is 
variable  to  permit  the  user  to  satisfy  specific  needs. 

In  addition,  an  analytical  solution  calculation  for  constant 
thermophysical  properties  has  been  introduced  into  the  program  to  check 
the  numerical  results.  Finally,  the  user  can  obtain  the  temperature 
distribution  as  a  function  of  time  in  tabular  and  graphical  forms. 


4.0  RESULTS 


4.1  Analytical  Solutions 

Schriempf  (Ref.  10)  gives  analytical  solutions  to  the  heat 
diffusion  equation  in  the  one-dimensional  situation.  He  assumes  the 
laser  beam  is  constant  and  uniform,  the  material  parameters  are  tem¬ 
perature  independent,  the  solid  is  uniform  and  isotropic.  Furthermore, 
there  are  no  internal  heat  sources  or  sinks,  and  no  phase  change  in  the 
materials  is  considered. 

The  solution  for  the  semi-infinite  solid 

t  /•_  «.'i  _  2.  e.  $.  COS  0  i  ■  v  1  .  f. 

T  (z,t)  =  - 1 -  r  0t  .  ierfc 

K  y 

and 

T  (0,t)  *  2.. COS  Q  gt[ 

K  fir 


hvH  ,m 

[22a] 


Similarly,  the  solution  for  a  slab  of  finite  thickness  is 
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T(z,t)  -  E^jcos  ej 


■2  -  L2 


.L.  (out  (L-z) 

— U2  L  ^ 

%  2  .  e-^2t/i2  .(cos  SlXtsiijl  | 


[23] 


and 


T(0,t)  = 

I 


£.<f>.  (COS  9] 


.  L.  fa.t  T  1  2  1  -an2ir2t/L2l ) 

— r*[  r^7‘  Jl 


[23a] 


where  L  is  the  slab  thickness,  and  the  condition  for  an  insulated  back 
surface  is 


[24] 


Equation  23  has  been  included  into  the  program  appearing  in 
Appendix  A,  such  that  the  reliability  of  the  numerical  results,  for 
similar  conditions,  can  be  assessed. 

4.2  Numerical  Solutions 


The  numerical  results  obtained  in  calculating  the  function  U  are 
identical  to  those  of  the  temperature  T,  if  we  assume  that  the  thermo¬ 
physical  properties  are  temperature  independent.  Otherwise,  we  have  to 
solve  the  following  equation 


U(T) 


1_ 

K, 


O 


K(T')  dT' 


[25] 


1 1,\'< .  LASS  11-11.:* 
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to  find  the  temperature  distribution  through  the  slab  considered.  For 
metals  studied  here,  the  thermal  conductivity  of  the  solid  phase  is 

K(T')  -  K.  (1  ♦  ef)  [26] 

Substituting  eq.  26  into  eq.  25  gives 

U  *  T  ♦  —  T2  [27] 

2 

and  finally, 


Furthermore,  we  assume  that  the  thermal  diffusivity  is  as 

follows : 

a(T.)  .  *1?')  .  Ml  *  BT-)  [29] 

C.  C  (T'J  P.C  (1  <■  tT1) 

P  P« 

where  K0  is  the  thermal  conductivity,  Po,  the  density  (assumed  constant 
in  the  temperature  range  considered)  and  Cp  ,  the  specific  heat  at  the 
initial  temperature  (i.e.  room  temperature  in  our  case). 


4.2.1  Numerical  vs  Analytical  Solutions 


This  section  compares  the  numerical  and  the  analytical  solutions 
for  similar  conditions;  it  also  gives  an  assessment  of  the  accuracy 
and  reliability  of  our  numerical  model. 
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Figure  1  shows  the  analytical  solution  (alternating  lines)  and 
the  numerical  results  (continuous  lines),  with  temperature  independent 
thermophysical  properties,  for  a  6-node  case  (i.e.  Az  =  0.2  cm)  along 
the  2  axis.  Each  curve  represents  the  temperature  distribution  at  5 
time  step  increments  (i.e.  5At) .  Both  models  predict  the  same  front 
surface  temperature,  within  1%,  after  only  25  time  steps  (25  At)  or 
about  0.2  s  in  the  present  case.  Furthermore,  if  we  calculate,  for  a 
specific  time,  the  area  under  each  curve  we  obtain  the  same  result,  which 
demonstrates  that  we  satisfy  the  energy  conservation  principle. 

The  same  results  appear  in  Fig.  2,  but  for  a  21-node  case  (i.e. 

Az  *  0.05  cm)  along  the  z  axis.  Each  curve  corresponds  to  16  time  steps. 
Except  for  the  first  few  curves,  it  is  very  difficult  to  differentiate 
between  the  results.  As  a  matter  of  fact,  we  have,  within  0.5%,  simi¬ 
lar  results  on  the  front  surface  after  only  16  time  steps  or  about  0.006 
s  in  this  case.  However,  this  accuracy  was  obtained  at  the  expense  of 

increasing  the  computation  time  by  a  factor  of  64.  This  is  because  Az 

2 

is  4  times  smaller  and  At  is  16  times  smaller  since  At  *  (Az)  . 

This  demonstrates  that  the  numerical  solution  can  give,  with 
some  compromise,  results  comparable  to  those  of  the  analytical  one.  In 
the  next  section,  the  numerical  solution  is  applied  to  cases  where  no 
analytical  solutions  exist,  but  where  there  are  a  few  approximate  ones. 
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4.2.2  Teaperature -Dependent  Themophysical  Properties 

The  temperature  dependence  of  themophysical  properties  is  an 
important  factor  in  dealing  with  laser-aatter  interactions.  Figure  3 
shows  some  results  for  copper  with  6  nodes  along  the  z  axis.  We  have 
assumed  that  the  themophysical  properties  vary  linearly  with  tempera¬ 
ture.  The  continuous  curves  represent  the  teaperature-dependent 
results  while  the  alternating  curves  are  the  results  for  the  constant 
room- temperature  case.  The  time  interval  between  curves  is  0.25  s. 

Both  temperature  distributions  are  similar  at  0.25  s  but  quite  differ¬ 
ent  after  4.5  s.  We  calculate  an  error  of  10%  or  about  0.5  s  of  the 
total  time  to  reach  melting  on  the  front  surface  between  the  two  models. 
In  other  words,  it  takes  0.5  s  longer  to  reach  melting  on  the  front 
surface,  for  the  specified  conditions,  if  you  consider  temperature- 
dependent  properties  as  in  the  physical  world. 

Figure  4  shows  the  absolute  error  in  temperature  for  the  prece¬ 
ding  conditions.  The  error  increases  from  the  front  to  the  back 
surface.  These  results  indicate  a  decrease  in  the  thermal  diffusivity 
of  copper  with  temperature  as  expected  (Ref.  11).  Figure  5  presents 
the  relative  error  in  percentage  for  the  same  conditions. 

Figure  6  illustrates  the  results  obtained  for  stainless  steel 
#  304.  This  case  demonstrates  the  drastic  influence  of  the  temperature - 
dependent  themophysical  properties.  Under  the  present  conditions, 
there  is  a  twofold  increase  in  the  time  required  to  reach  melting  on 
the  front  surface  between  the  two  models.  In  the  variable  model, 
melting  on  the  front  surface  can  be  reached  in  3  s,  whereas  in  the 
constant  model,  this  can  be  achieved  within  1.5  s. 
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FIGURE  3  -  Temperature  distribution  within  a  1  cm-thick  Cu  target  (6  2 

nodes  along  z  axis)  for  an  incident  intensity  of  20,000  W/cm 
and  a  coupling  coefficient  of  0.04.  The  time  interval 
between  each  curve  is  0.25  s  for  the  temperature-dependent 
case  (continuous  lines)  and  for  the  constant  room-temperature 
case  (alternating  lines) . 
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FIGURE  S  -  Percentage  of  error  in  temperature  distribution  between  the 
two  models  of  Fig.  3. 
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FIGURE  6  -  Temperature  distribution  within  a  1  cm-thick  stainless  steel 
#304  target  (21  nodes  along  z  axis)  for  an  incident  intensity 
of  20,000  W/cm^  and  a  coupling  coefficient  of  0.04.  The  time 
interval  between  each  curve  is  0.25  s  for  the  temperature- 
dependent  case  (continous  lines)  and  for  the  constant  room- 
temperature  case  (alternating  lines). 
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Furthermore,  the  curves  cross  one  another  and  we  find  a  higher  value  of 
the  temperature  on  the  back  surface  in  the  variable  case.  Figure  7 
shows  the  absolute  error  in  temperature  along  the  z  axis  for  the  same 
time  interval . 


4.2.3  Averaging  of  the  Thermophysical  Properties 


We  define  the  average  value  of  the  thermophysical  properties  in 
the  following  way: 


K 


av. 


K(T.)  +  K(Tm) 
2 


and 


C 

pav. 


C  (T.)  +  C  (T) 
p'-  i/  pv  -’m 


[30] 


[31] 


where  is  the  initial  temperature  and  the  melting  temperature.  If 
we  use  these  values  in  the  constant  properties  model  instead  of  those  of 
the  room  temperature,  we  obtain  the  results  shown  in  Fig.  8  for  copper 
and  those  in  Fig.  9  for  stainless  steel  #304.  We  expect  similar  results 
because  the  thermophysical  properties  of  these  two  materials  vary  linear¬ 
ly  with  temperature  so  that  these  properties  are  overestimated  at  low 
temperature  and  underestimated  at  high  temperature  by  the  average  pro¬ 
perties  values.  However,  although  the  discrepancy  in  time  to  reach 
melting  on  the  front  surface  is  removed,  the  temperature  distribution 
through  the  slab,  specially  in  stainless  steel  #304  is  quite  different. 
The  distribution  behaves  as  if  the  coupling  coefficient  was  varying  with 
temperature  and,  depending  on  the  position  selected  along  the  z  axis, 
it  seems  to  increase  or  to  decrease  with  time. 
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FIGURE  7  -  Absolute  error  in  temperature  distribution  for  conditions  of 
Fig.  6. 
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FIGURE  8  -  Temperature  distribution  within  a  1  cm-thick  Cu  target  (6  2 

nodes  along  z  axis)  for  an  incident  intensity  of  20,000  W/cm 
and  a  coupling  coefficient  of  0.04.  The  time  interval  between 
each  curve  is  0.25  s  for  the  temperature -dependent  case 
(continuous  lines)  and  the  average-value  case  (alternating 
lines) . 
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FIGURE  9  -  Temperature  distribution  within  a  1  cm-thick  stainless  steel 
#304  target  (21  nodes  along  z  axis)  for  an  incident  intensity 
of  20,000  W/cnr  and  a  coupling  coefficient  of  0.04.  The  time 
interval  between  each  curve  is  0.25  s  for  the  temperature- 
dependent  case  (continuous  lines)  and  the  average-value  case 
(alternating  lines) . 
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FIGURE  10  -  Temperature  distribution  within  a  0.1  cm-thick  A1  target 

(21  nodes  along  z  axis)  for  an  incident  intensity  of  2  x  10° 
W/cm2  and  a  coupling  coefficient  of  0.02.  The  time  interval 
between  each  curve  is  134  x  10’  s  for  the  numerical  case 
(continuous  lines)  and  the  analytical  average  value  case 
(alternating  lines) . 
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Finally,  Fig.  10  presents  some  results  for  aluminum,  using  the 

average  properties  values  as  defined  previously.  We  used  21  nodes  for 

this  0.1  cm-thick  slab  because  we  were  interested  in  a  high-input 

6  2 

intensity  ($  =  2  x  10  W/cm  )  for  a  short  time  duration,  as  would  be 
encountered  with  pulsed  lasers.  The  time  to  reach  melting  on  the  front 
surface  is  within  1%  of  the  one  obtained  by  Gautier  (Ref.  12). 


S.O  DISCUSSION 

Our  aim  of  developing  a  reliable  numerical  procedure  for  the  non¬ 
linear  one-dimensional  heat  conduction  equation  has  been  achieved.  The 
boundary  conditions  are  calculated  through  a  numerical  procedure  differ¬ 
ent  from  the  one  used  by  Hanley  (Ref.  1) ;  therefore,  we  satisfy  the 
energy  conservation  law  whereas  his  scheme  does  not.  Furthermore,  our 
scheme  is  simpler,  easier  to  understand  and  to  work  with,  and  through 
proper  selection  of  the  "stability  constant",  converges  more  rapidly  to 
the  true  solution  for  an  optimal  number  of  steps.  The  time- independent 
or  continuous-intensity  laser  beam  has  been  considered  in  the  present 
work. 


Our  present  solution  is  applicable  to  cases  where  lateral  conduc¬ 
tion  of  heat  is  negligible.  These  situations  correspond  to  short  inter¬ 
action  time  consideration  or  to  experiments  in  which  the  laser  beam 
dimensions  are  similar  or  larger  than  those  of  the  target.  However, 
extension  of  our  numerical  model  to  3  dimensions  is  easy  and  straight¬ 
forward  . 

Similarly,  the  model  does  not  include  any  explicit  terms  for 
radiative  and  convective  heat  losses  but  these  could  be  considered  by 
modifying  the  appropriate  boundary  conditions. 
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6.0  CONCLUSIONS 

A  numerical  method  of  solving  the  nonlinear  partial  differential 
equations  of  heat  conduction  with  a  computer  has  been  developed.  The 
reliability  of  the  numerical  results  has  been  demonstrated  by  compari¬ 
son  with  the  analytical  solution  in  one  dimension.  By  proper  selection 
of  the  "stability  constant"  the  numerical  procedure  converges  more 
rapidly  to  the  results  of  the  analytical  solution  by  minimizing  the 
truncation  error.  Agreement  is  within  1%  after  only  25  time  steps 
(25  At). 

We  have  shown  the  importance  of  taking  into  account  the  effect 
of  temperature -dependent  thermophysical  properties  when  dealing  with 
laser-matter  interactions.  For  example,  in  the  case  of  stainless  steel 
#304,  the  calculated  time  required  to  reach  melting  on  the  front  sur¬ 
face  in  the  case  of  variable  properties  is  twice  the  time  necessary 
when  room -temperature  values  are  used.  Although  this  discrepancy 
disappears  when  average-temperature  values  are  used,  the  calculated 
temperature  distributions  through  the  material  are  significantly  dif¬ 
ferent  . 
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APPENDIX  A 

A. 1  Program  Description  and  Listing 

The  program  written  in  Fortran  IV  may  be  used  under  the  CP-V 
Version  E01  operating  system  on  DREV  Xerox  560  MP  computer.  The  listing 
shows  the  materials,  with  their  corresponding  thermophysical  properties, 
available  to  the  user.  Having  selected  the  material,  the  user  inputs 
the  thickness  of  his  target,  the  number  of  nodes  through  it,  the  ambient 
temperature,  the  incident  flux,  its  angle  of  incidence  from  normal,  the 
coupling  coefficient,  the  total  time  of  run,  the  print  out  interval  and 
the  type  of  boundary  conditions.  From  there  on,  the  program  starts  its 
calculation  by  establishing  the  initial  conditions  throughout  the  slab. 
Then,  it  determines  the  temperature  for  the  next  time  increment  and 
checks  if  it  has  to  print  out  those  results.  In  the  event  it  has,  we 
obtain  tables  of  values  for  the  front  and  back  surface  temperature  and 
its  distribution  through  the  slab.  Furthermore,  it  calculates,  for  this 
time,  the  corresponding  analytical  solution  and  outputs  it.  Then,  the 
program  checks  if  the  total  run  time  is  exhausted  and  if  not,  it  goes 
back  through  the  loop  and  calculates  the  temperature  distribution  for 
the  next  time  step.  Once  the  alloted  time  is  passed  the  program  calls 
subroutines  to  create  APL  files  with  the  complete  set  of  data.  These 
files  are  used  to  produce  graphical  representations  similar  to  the  ones 
appearing  in  this  report. 
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1.000 
2.000 
3 . 000 
5  .  500 
6 . 000 
7  .000 
8 . 000 


o JOB  66043 .JPBMORENCY ,7 .  ( TERMINAL  JOB  BFALACOC A ) . 

•  LIMIT  (TIME ,  20)  ,  (CORE  t^8) 

o ASSIGN  F : 1 05 , (FILE , HTINPUT1 ) , (IN) 

•  ASSIGN  M: SI , (FILE , PALACOC  A 1 ) , (IN) 
o FORTRAN  LS,SI ,G0 ,BC 

* LOAD  (GO) , (LMN ,FALACOCALM) ,(PERM) , (EF , ( APLFNS , LPR) ) 
o RUN 
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ONE-DIMENSION  HEAT  CONDUCTION  LISTING 


1.  DIMENSION  U010:21*2i>T0<2i,a),T3f**(81il50),T0<2i  l#TM!2l#100> 

2*  DOUBLE  PRECISION  PS 

3*  I  FORMAT! •  1'  > 

*•  10  FORMAT! '  1  •••ALUMINUM#  2  • * • COPPER*  3  • • •MOLYBDENUM* ’ I 

5.  20  FORMAT!'  *  ...IRON#  5  •••NICKEL#  6  ...TITANIUM#  7  • • . ALl 202* 1 ' 1 

6 *  25  FORMAT! •  8  « . ‘STAINLESS  STEEL ( 30* )• ) 

7.  30  FORMAT!'  INPOT  8  FOR  TARBET  MATERIAL' ) 

8*  <.0  FORMAT (  I* ) 

s •  so  format (•  material  thickness  icm)'> 

10*  51  FORMAT! 'NUMBER  OF  NODES  ALONG  Z-AXIS') 

11*  55  FORMAT! 2I5#F20. 51 

12«  60  FORMAT ( F20 .5  I 

13*  90  FORMAT!'  AMBIENT  TEMPERATjRE  (O') 

I*.  ioo  format t •  incident  intensity  ihatts/cm*cm)  »> 

15*  120  FORMAT ( '  ANOlE  OF  INCIDENCE  IDEG)  FROM  NORMAL' I 

16*  123  FORMAT!'  PREmElT  ABSORPTION  COEFFICIENT') 

17.  130  FORMAT!'  TIME  OF  RUN  (SEC)') 

18*  135  FORMAT!'  PRINT  OUT  INTERVAL  (SEC)') 

19*  150  FORMAT!'  1  •••BACK  SURFACE  INSJLATED  ') 

20*  lbO  FORMAT!'  2  •••BACK  SURFACE  HEAT  SInkEO  TO  AMBIEnT  ') 

21*  185  FORMAT! 'FRONT  SURFACE  STARTS  TO  MELT  AT  POSITION  ll> 

22*  1  AND  AT  TIME  T"  ' #F10»5# '  SECONDS') 

23 •  190  FORMAT! ' 1TIME  IN  SECONDS '#  PU  .7# // )  • 

2%.  191  FORMAT! 'OFRONT  SURFACE' ) 

25*  192  FORMAT ( rF15  <8 ) 

26<  193  FORMAT! 'O0ACK  SURFACE') 

27.  19*  FORMAT! 'OTEMPCRATURE  ALONG  Z-AXlSM 

28*  195  FORmaTC  RUN  COMPLETED') 

29»  196  FORMAT! 'OanAL^TICAL  SOLUTION  ALONG  Z-AxIS») 

30 •  WRITE  ( 108# 1 ) 

31  •  C  SELECT .  MATERIAL  code 

32 •  WRITE  Il08#10> 

33*  WRITE  ( 108#  20  I 

3**  WRITE  ( 108# 25  1 

35*  WRITE  I  1 08#  30 1 

36*  READ  ( 108#  *0 1  111 

37.  OUTPUT  Ill 

3a •  GO  TO  »200#21Q#220#230#2*0#2*1#2*2.2*5)# III 

39 •  200  CONTINUE 

*0«  c  Ill«l  •••ALUMINUM 

*!•  A0»660. 

*2*  F0« • 25 

*3*  F2». 57361 

**•  A3"2 • 7 
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45* 

*4-94. 

46« 

Fl-O. 

47. 

F3-0. 

*8* 

F4-0. 

49. 

F5-30  • 

go. 

F9-2.43 

bl» 

30  TO  300 

52. 

210 

CONTINUE 

b3* 

C 

111*2  ...COPP£« 

54. 

*0*1083* 

65* 

F0-.091 

66* 

F2-.956 

67. 

*3-*. 89 

68  * 

*4*42* 

69* 

F1-2.456E-5 

60* 

P3**1.567£*4 

fci  • 

F4*0« 

62* 

F5-0.0 

63* 

F9-8.217 

64. 

30  TO  300 

65* 

220 

CONTINUE 

66* 

C 

111*3  ...mOlYBD 

67. 

*0-2610* 

68* 

F0-. 06162 

69* 

F2-.3460 

70* 

*3-10.2 

71. 

*4-131 . 

72* 

F1-2.2E-5 

73. 

F3  —  3.46E-5 

74. 

F4-0. 

75* 

F5-10. 

76* 

F9-8.1 

77. 

30  TO  300 

78* 

230 

CONTINUE 

79* 

C 

111-4  . • • I HON 

80. 

*0-1535* 

el  • 

FO- . 1060 

e2  • 

F2-.1080 

83* 

*3-7.88 

84. 

A4-65* 

85. 

F1-9.6E-5 

86* 

F3—  1.08E-6 

87. 

F*-0* 

bb< 

F5-18. 
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89* 

F9-6.88 

90. 

30  TO  300 

91  • 

2-0 

continue 

92* 

C 

111-5  ...NICKEL 

93* 

AO-1453* 

94  • 

F0-* 1095 

95. 

F2-.1485 

56* 

A3-8.7S 

97. 

A4-73. 

98* 

Fl-5«  49£-b 

99* 

F3--4 . 56E-5 

1  0&  • 

F4-0« 

Ifcl* 

F5-0  • 

1C2  • 

F9-7.9 

103* 

ao  To  300 

104. 

241 

CONTINUE 

105  • 

C 

111-6  ...TITANIUM 

106* 

AO-1690* 

1C7. 

F0-.139 

108* 

F2-.0372 

109" 

A3-4.54 

110* 

A4-103.9 

Ill* 

Fl-O. 

112* 

F3--4 • E“6 

113- 

F4-0  . 

114. 

F5-21 • 

115* 

FS-4.09 

1 1  8  * 

80  TO  300 

117. 

242 

CONTINUE 

118" 

C 

111-7  •  *  *  AL ( 2024  > 

119* 

A0-630 . 

120. 

F0-.215 

121* 

F2-.33- 

122* 

A3-2.77 

123* 

A4— 95 • 6 

124. 

F1-7.7E-5 

125* 

F3-9.0E-4 

126* 

F4-0 . 

127. 

F5-0. 

128. 

F9-2.43 

129* 

ao  TO  300 

130* 

245 

continue 

131. 

C 

111-8  . . >aCIER(304I 

132. 

AO-1450* 

ai 
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133*  F0«. 11712 

13*»»  F2*.  03615 

135*  a3»7*9 

136*  A4«65« 

137*  Fl»3.7<.E-5 

138*  F3-3.32E-5 

139*  F 4*0  « 

140*  F5«20 • 

141*  F9«7. 

Iw2.  300  CONTINUE 

1<*3»  C  INPUT . MATERIAL  THICKNESS  (CM) 

144.  WRITE  ( 108/ 50 ) 

1  *.5  •  READ  (  105/60  !  A6 

146*  OUTPuT  A6 

i*»7.  c  input . number  of  nodes  alons  e-axir 

l48»  WRITE! 108/ 51 ) 

lw9*  READ ( 105/ 60  )  KK 

ISO •  OUTPuT  KK 

151 «  NK«*KK-1 

162*  C  INPuT. ....  AMBIENT  TEMPERATURE  (Cl 

153.  WRITE  ( 108/ 90 ) 

16**  •  READ  (  105/60)  A9 

155 •  OUTPuT  A9 

156 •  c  input .  incident  intensity  i watts/cm*cm 

157.  WRITE  (108/100) 

168 •  READ  (105/60)  BO 

159*  OUTPUT  80 

lbO  •  B0*B0/4 • 18a 

1 6 1 •  C  INPUT . angle  of  incidence  (DEQ)  from  normal 

lb2 •  WRITE  (108/120) 

163.  READ  (105/60)  62 

164.  OUTPUT  82 

165*  C  INPuT .  PREMELT  ABSORPTION  COEFFICIENT 

166>  WRITE  (108/123) 

167.  READ  (105/60)  A5 

168*  OUTPUT  A5 

169.  C  INPUT .  TIME  OF  RUN  (SEC) 

170*  WRITE  1108/130) 

171.  READ  (105/60)  B3 

172.  OUTPUT  83 

173*  C  INPUT. ....  PRINT  OUT  INTERVAL  (SEC) 

174.  WRITE  (108/135) 

175.  READ  (105/60)  DO 

176.  OUTPUT  DO 
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177.  C9-0. 

178*  C  SELECT . BOUNDARY  CONDITIONS!  195 

179*  WRITE  I  108.  150 ) 

lbO •  WRITE  (108/160) 

1 B 1  *  . . .  •>  BACK  SURFACE  INSULATED 

ib2*  . . .  •)  back  surface  heat-sinked  to  ambient 

183.  READ  (105/40)  185 

184.  OUTPUT  IB5 

185*  C  FRONT  SURFACE. •••• INCIDENT  FLJX 

186*  P-80 

1 87 •  01-00 

1 68 •  INDI-O. 

1B9<  C  INITIAL  conditions 

190*  DO  344  K-l/KK 

191*  TO ( K/ 1 1 ■ A9 

192*  IF ( F3 .£0 * 0  *  I  30  TO  345 

193*  uO  (  </  1  )-A9-F5-F3»F54A9/'F2*.5*F3*F5*F5/F24.5AF3*A:»*A9/F2 

194.  QO  TO  344 

195*  345  UO ( K / 1 ) -A9 

196*  344  CONTINUE 

197.  CTE-1./6. 

198*  C2-CTE*  <F04F14A9)»A3aa6*A6/(F24F3*(a9»F5> l/NK/NK 

199*  C0-C2*NK*NK/A6/A6 

200 •  C4-2>4C2*NK4nK/A6/A6 

201.  THET-3. 141592*82/180. 

2u2 •  33-C0S(ThET)4A5*A6*2*/NK 

203*  500  IF(T0(1/1)«LE.A0)  00  TO  505 

204.  WRITE! 108/ 185I/C9 

2C5 •  QO  TO  1000 

206-  505  CONTINUE 

2C7.  F8-A3 

208«  00  560  K-l/KK 

209*  U0(0/1)-U0(2/1)4Q3»P/'F2 

210.  QO  TO  ( 530/535 ) / IB5 

211*  530  UO(KK4i/l)-uO(KK-l/l> 

212*  QO  TO  550 

213*  535  U0(KK41/1)-2.4U0(KK/1)-U0(KK»1/1) 

214.  550  CONTINUE 

215*  Al-<F2*P34(T0(K/i>-P5> ) /FB/ ( F0*F1 »T0 ( <, l  )  ) 

216*  El-Al*C04(U0(K4l/li»U0tK-l/i) ) 

217.  U0!K/2)-El4( i.-Al*C4)6U0(K/i ) 

218*  IF ( F3 . EQ • 0 • 0 )  QO  TO  555 

219*  BETA-F3/F2 

220.  TO  ( K/  2  )  -  (  <  (  ( l«42«*3ETA*U0(</2>  >  »*  .  5  > /BETA  )  4-F5-1  * -'BETA  ) 
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221*  00  TO  560 

222 «  555  TO  I </ 2  I -U0 ( to  2 ) 

223*  560  CONTINUE 

224.  C9-C9+C2 

225*  580  IF i C9  «lT  *  01 )  00  TO  900 

226*  01-01400 

227.  INDI-InDI*! 

228.  00  600  K-l  / KK 

229*  TGRa ( K / InOI I -TO  I K/ 1 ) 

230.  600  CONTINUE 

231*  TDIM-C9-C2 

232-  WRITE! 108/ 190)  TDI* 

233.  WRITE ( 108* 191 1 

234 •  WRITE! 108/192)  (T0I1/1>> 

235.  WRITE! 108/193) 

236*  WRITE! 108/192)  (T0(KK/1>> 

237.  WRITE ! 108/ 194  I 

238>  WRITE ( 108/ 192 )  < TO ! 1/ 1  I / L-l/ KK ) 

239*  . . ANALYTICAL  SOLUTION 

240*  LL-lNOI 

24 1 •  PI-3.14159265 

242*  0IFF-IF24F3*! A0-F5J/2. ) /A3/ l FOtPI * < A0-F5  )  /2. ) 

243*  700  Z-0.0 

244.  J-1 

245*  725  wl-EXPl-0lFF*PI*Pl*TDIM/A6/A6) 

246*  N— 1 

247*  C-0.0 

248*  750  W-l-l )-4N/N/N*EXP(-0IFFan*N*PI4PI*TDI'1/A6/a6) 

249*  W2-W*C0S(N*PI*(A6-Z)/A6> 

250*  Vl-ABSIW/Wl) 

251*  IF!vI.lE.O. 000002)  30  TO  775 

252.  C-C4W2 

263*  N-N+l 

254.  QO  TO  750 

255*  775  C-C+W2 

256*  TO! J)-a5*P»TDIM4DIFF/!F24F3*( A0-F5I/2. ) /A64A5*p*A6/ I F2*F3* ( 

257.  1 A0-F5 )/2* )*! ( 3* ( A6-Z > **2- a6*A6 I /6 • / A6/A6-2 . *C/P I /P I )*A9 

258*  Z-Z4A6/20. 

259*  J-J+l 

260*  IF(Z.ar«A6)  GO  TO  800 

261*  GO  TO  725 

262*  aOO  WRITE! 108/196) 

263>  WRITE  1 108/ 192 )  I  TO < U ) / U-l / 21 ) 

264.  DO  850  1-1/21 


1 

i 

i 

! 

1 

\ 

1 


i 


\ 

I 
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EtoS  « 
266* 
Cfe?* 

EoS  • 
ce.9  • 
2  70  • 
£  7 1  • 
w72  • 
£  73  • 
27**  • 
£75  • 
E7t>  • 
£77. 
27S  • 
279  • 


T,-(  I,Ll>»T0(  I  1 
e50  CONTINUE 

. . END  OF  ANALTTJCAL  SOlJTIOn 

SuO  IFIC9*aT«B3*0**C9»EQ»B3>  30  TO  950 
DO  9 1 0  K«liKK 
uO ( K> 1 ) “UO ( K<  2 ) 

TQ(  n> 1 |»T0(K<2) 

910  CONTINUE 
GO  To  500 

350  WRITE ( lOdj 195 ) 

1000  CONTINUE 

CalL  GwAPM(TriJlNDI*21«l  > 

CALL  QSAPhITGRa* INDIa«j 5) 
call  exit 

END 
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2 
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4 

5 

6 
7 
b 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 


SUBROUTINE  GRAPH [TEMP* INDI*KK*IAN> 
DIMENSION  TEMPI«K*INDI ) 

DIMENSION  ITYPE(2)*ISIZE(3) 

I  SIZE  < 1 ) *2 
ISIZE ( 2 ) *21 
ISIZE ( 3 ) "INDI 
ITYPE<  U-3 
ITYPE<g>*4 

CALL  FTIE< 1* ' JCAOUT' > 

CALL  FREPLACE ( 1* I  AN* TEMP* ISIZE* I  TYPE  ) 
ITYPEt 1 )"2 
ITYPE(2)-2 
ISIZE ( 1 ) "0 

CALL  FREPLACE Jl* 100* <<* ISIZE* I  TYPE) 
CALL  FREPLACE t 1*101* I VDI* ISIZE* I TYPE) 
CALL  FUNTIEI1) 

RETURN 

END 
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A. 2  Input  and  Output  Examples 

The  program  is  run  from  a  terminal  in  two  possible  modes.  In 
the  ON-LINE  interactive  mode  the  program  asks  for  the  specific  input 
file  considered  and  then  starts  the  calculation.  The  output  is  dis¬ 
played  on  the  terminal  used  and  is  similar  to  that  obtained  in  the 
BATCH  mode.  In  the  BATCH  mode,  a  job  file,  which  contains  a  set  of 
controls,  the  program  and  data  input  and  output  files,  are  run.  The 
output  is  printed  on  the  computer  line  printer  as  shown  in  the  typical 
example  shown  hereafter. 
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INPUT  DATA  FILE 


1  •••ALUMINUM.#  2  •••COPPER*  3  •  • 'MOLYBDENUM# 

*  ...IRON#  5  ...NICKEL#  6  ...TITANIUM#  7  ...AH202V) 
g  •••STAINLESS  ST£ELOO%I 
Input  *  for  Target  material 
hi  •  a 

material  Thickness  (Cm) 
as  •  i«ouooo 

NUMBER  OP  NODES  ALONG  Z-AXlS 
KK  ■  21 

AMBIENT  TEMPERATURE  Id 
A9  *  20*0000 

INCIDENT  INTENSITY  (WATTS/CMACM) 

BO  ■  20000*0 

ANGLE  OF  INCIDENCE  (DEG)  FROM  NORMAL 
62  •  *000000 

PREhELT  ABSORPTION  COEFFICIENT 
as  -  A.OOOOOOE-02 
TIME  OF  Run  (SEC) 

63  »  3*00000 

PRINT  OUT  INTERVAL  (SEC) 

DO  ■  *250000 

1  ...back  surface  insulated 

2  ...BACK  surface  heat  bInKED  TO  AMBIENT 
IBS  ■  1 
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OUTPUT  DATA  FILE 


time  in  seconds  •a»6S*s3 


front  Surface 
5l?*17797a52 

SACK  SURFACE 

£C*OOCOOuOC 


TEMPERATURE  AL 
Sl7.17797B52 

£7 • 67676787 

£0  *00122070 
EC  ■ OCCOCwOG 


Z-AXlg 

3a9  •  e0bISa3G 

22*62307688 

20 • COOGUOOO 
20 * OOOGuOGO 


217*79727305 

20*77685577 

20.00000000 

20*00000000 


125*55175781 

20.20263672 

20.03000000 


analytical  solution  AlOnS  Z-aXIS 


*30*66012695 
32*62736965 
£C  *02*26177 
15*99971008 


290*82963398 

25**S023967 

20*00599670 

19*99990875 


188*17067026 

22*19616699 

20*00085779 

20*00009155 


117*67336257 

20*81092837 

19*99990875 


70.03808597 

£0.0*687500 

£0*30000000 


73*2001 1902 
20*275*3670 
19*99990875 


71.05175781 

20*00952178 

20*00000000 


76*90753177 

20*08557129 

19*99990875 


time  IN  SECONOS  *7936965 


front  surface 
693*62307688 

SACK  SURFACE 
co 'OuOOOuOO 


TEMPERATURE  ALI 
693*6230*688 
6**3371 5b20 
£0*6*916992 
£0* OuOOOuOO 


Z-AXlS 

537*38012695 
53*30810577 
20.2*731775 
20 • OOOuUOQO 


399*82753906 

37*02905273 

20*0|837b91 

20*00000000 


283.76118167 

28*20971680 

20.02905273 


1 92  *53772383 
£3*73713086 
20 ■ 00952178 


126*16528320 

21*60327178 

20.00219727 


analytical  SOLUTION  alOnS  Z«aXI8 


601 *07006836 
8**  *  19859317 
21*93779*02 

20*0127*109 


*55.33862305 

59***902039 

20*93702698 

20.00*77578 


336*76177761 

73*37881592 

20*73501282 

20*00303650 


273.76356201 

33.29827881 

20.19317575 


172*60578918 

£7*2b207397 

20*06219910 


120.73976135 

23*83071392 

20*03311157 


XHISTAGSI5E1 

Gu-  1 
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TImE  IN  SECCNOi  .7*05*47 


FfiflNT  SwRFACE 

82u*00000u00 

back  SURFACE 
20  *01562500 

TEmPlRATuRE  ALONG  Z-AXlS 

820 .00000-00  672.2712*023  536.356398** 

159*  SS891602  110.8*252930  76.50537109 

26.6230*6 8b  22.86*01367  21. *0307617 

20*05517578  20.02*9023*  20.01562500 

analytical  SOLUTION  alOnO  Z-AXI8 

731.66186523  583*39*0*297  *57*50317383 

1*6*86979675  109.708*1980  80.96261597 

jo. 10*62552  26.03887939  23.51322937 

20*3215179*  20.19039917  20*15176392 


TI«E  IN  SECONOS  .9981252 


FPnM  Surface 

925*62231**5 

back  surface 
2C  >2585**92 

TEaRERATuRE  ALONG  Z-AXlS 

925.62231**5  78*.702l*a**  652.3210**92 

2*7.17163086  162*8*936523  132.98091797 

38*9*7021*8  31.23779297  26>*7*60938 

20*56298828  20.3271*8**  20.2565**92 

ANALTTICAL  SOLUTION  along  Z'AXIS 

8*8.210**522  696.37353516  965.90307617 

218*6***85*7  168.07295227  128**738*6** 

*£.0*951*77  37.39862061  31.*02*0*79 

21*89921570  21.3*771729  21.17568970 


*15.2*316*06  311*391357*2  226.2580566* 

53.70288086  39*3000*883  30.62207031 

20.66162109  20  *29907227  20*13085938 


352.9**3359*  268*06398*38  200.8*373*7* 

60. *18*1125  *6*13006592  36*5631805* 

21.969*8669  21*0977*780  20.59326172 


530.60986328  *21*55859375  326.7661132s 

95.79077168  69*2*682617  51*013*2773 

23.62622070  21*9/729*92  21.05395508 


65*. 12359516  359*96093750  282.01391602 

98.0625*578  75.16*55078  58*2665710* 

27.33912659  2**65*55627  22.93768311 
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TI(-E  IN  SECONDS  1  •  2*«95*1 


FRONT  SURFACE 
loll*  3*6*3555 


sack  surface 
21*30*19922 


TEpPLRATuKE  alOnq  Z-AXIS 
1011 *3<*6*3b55  875*65283203 

32t  • 59087695 


56*89379883 
cL • £7929bs8 


257*606889*5 
*5 • 90*29688 
21*53271*6* 


7*6*6*1357*2 
193*66650391 
36 **£2607*2 
21*30*19922 


625*933105*7 
1*5*12255859 
33*1 7968750 


513*12182617 

107*7*9023** 

26*1b79B82B 


ANALYTICAL  SOLUTION  ALOnO  Z-AXIS 


9*2.730*6875 
28* • *70*5898 
67*03033**7 
25*65603210 


791*87182617 

225**110107* 

53*53*17969 

2***0707397 


658**01855*7 

tFT'37393l8S 

63*58168030 

*4*01371765 


5*1**3129883 

138**97018*3 

36*38682556 


**1*37695313 

108*33567932 

31*31039*29 


TI*iE  IN  SECONDS  I**9l7s31 


FRONT  SURFACE 
1066*60*2*805 


SACK  SURFACE 
-  a2*  *012*51 17 


TEMPERATURE  ali 
1086*60*2*605 
.63*69165039 
67 • 6s380b59 
26*12*51172 


i  Z'AXlS 
955. 2607. 219 
325.27392378 
66*80*59375 
2**51928711 


829*37280273 
256*37207031 
51*73193359 
2* *012*5117 


710.21801758 

198.9155273* 

*1.12133789 


599*0920*102 

132*33009766 

33*8691*063 


ANALYTICAL  SOLUTION  AlOnQ  Z-AXIS 


1030*0686*766 

3*6*32910156 

92*12615967 

32*22*2*316 


878**5361326 
282  *  **6*8633 
73**1950376 
30.1*996191 


7*2.752**1*1 

227***bl2122 

59*7621*600 

29**7726**0 


622*59082031 

182.10*73633 

*9.1698608* 


517*36*99023 

1*5*20967102 

*1*38*86770 


IHIS  Pal-,  i-  ■  ;  •“ 

FHOM  GOP*  / 


lj'.i  *  r 

10  log 


*15.6*18*570 

79*93212891 

23*72265623 


356*00170896 

85*00825500 

27.851501*6 


*97  *  2268066** 
115*69091797 
29*09033203 


*26*26879883 

115*5831**0* 

35*8728*851 
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TIm£  IN  SECONDS  1*7*93*39 


FRONT  SJRFaCE 
1156*8686523* 

eiCN  SURFACE 
29*5*296875 

TEmRERATuRE  ALONG  Z“AXlS 
1156*8686523*  1029.38769531 

*81.09*97^70  396*23071289 

122*1799316*  9**2560S39* 

33*2563*766  30***55966* 


906 • *60**922 
321 **0991211 
73.075*39*5 
29*5*296875 


789*120605*7 

256*88818359 

57.36938*77 


ANALYTICAL  SOLUTION  alOnG  Z-AXIS 
1113*79S89o**  961*56835937  82**05737305  700.9670*102 

*12*63110352  3*0.91357*22  279*8*057617  228*39118958 

121*521255*9  98. *767*086  80*2616*2*6  66.13*15527 

*2 • *5678711  39**28*6680  38**3872070 


67a**2a95508 
202*56200195 
*6  *  08593750 


5? 1 • 8 1 660352 
105*52926636 
55  *  *5605*69 


575**3090820 

157.9199*531 

36*29370117 


*95*967773** 

150*23272705 

*7.70117186 


T  Im£  IN  SECONDS  1*9961729 


FRONT  SURFACE 
1216*06225586 

back  surface. 
38*1*111328 


TEMPERATURE  ALONG  Z-AXis 
2218*06225566  1093*79026320  973**3701172 

5*6*73120117  *61*0*028320  382*1799316* 

155*09350586  12**9565*297  98*02758789 

*3*65087691  39**9096680  38*1*111328 

ANALYTICAL  SOLUTION  ALOnQ  Z-aXIS 
1188*  *2285156  1035.72387695  896*81699*53 

*72*06201172  395*66967773  329*6721191* 

151*98117065  12**76715089  102*82881165 

55*511*7*61  51*5717773*  50*27758789 


857.8625*683 

312.55159375 

77.3*22*516 


771 .*587*023 
273.197021*8 
85. **120117 


7*7.9611816* 

232.32128906 

61*9*262695 


659*25537109 
22S • 3*90**80 
72*12719727 


6** • 6286621 1 
201.315*2969 
50*99*10156 


559*57358398 

195*23339318 

62*26931055 


\  '  O' 
,  .  'y' 


,<cv',Vv 


mmmm*  ~ 


UNCLASSIFIED 

46 


TjnE  IN  SECONDS  2*2*30019 


Fknf.T  SURFACE 
1?7*  **30*1992 

eACK  SUhFj,C£ 

5L • *0991211 


TEMPtRATyRE  ALI 
1?7***30*1592 
612*23535156 
190*568**766 
57*62*3691* 


Q  Z-Axls 
1152.99023*3* 
522 • 6679**75 
153.8857*219 
52.2379*8*3 


103*. 98266602 
**0.906**727 
126.6772*609 
50 >*0991211 


921.13525391 
367. **189*53 
101.22167969 


*12*21191*06 
302  >556396*8 
Bl • 7**1*063 


ANALYTICAL  SOLUT 
1 ?5a • 57**6289 
529 • *91*550* 
1B*«15600lb6 
71 .6125*bb3 


ION  AlOnU  Z*aX 
1105**8559570 
•*9.1219*203 
153.1679*72* 
66.80566*06 


965**20896** 

378.89770501 

127*82128906 

65*221191*1 


638.177*6023 

316*06787109 

107**9*62891 


723**2*56055 
265  *89*00391 
9 1 >6*18*570 


TI*E  IN  SECCNOS  2**898310 


FRnl.T  SURFACE 
1326*815*2969 


(JACK  SURFACE 
66-5*907227 


TEMPERATURE  AU 
1326*815*2969 
672.033**727 
233  *  732*2168 
75*66*7*609 


0  Z-AXI9 
1207.90591758 
581.25219727 
195.26977539 
68.85*2*605 


1092.0*3*5703 
*97 .*29**336 
158*3686523* 
66*5*907227 


979*657*2168 
*21 .073730*7 
128.5*5*1016 


872*00585937 

35**96229566 

109*23315*30 


ANALYTICAL  SOLUTION  ALOnQ  Z-aXIS 
1 32*  *  9bb2s1 25  1171*57031250  1030*52929688 

5*5  *096*3555  501.30*105*7  *27. *3017578 

217*68812561  183.3157806*  15**8918*570 

90*5*61*258  8**936923**  83*06276367 


901*69262695 
362.81009659 
131 .697*2188 


7o* • 7  7b56**5 
306*7*953125 
113*71*99*1* 


708*990*7*52 
2*6*3232*219 
67  *  *96093 79 


620.71313*77 

221**7293091 

79*80322266 


769*17065*30 

292*10693359 

87.8*987695 


679**0356**9 
258*621337*9 
130*  05961**5 
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TImE  In  SECONDS  2«7»739i7 


FHCinT  S-*FaCE 
1377.909*2363 
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physical  properties  variation  on  the  temperature  distribution  within  a  physical  properties  variation  on  the  temperature  distribution  within 
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